'''
Created on Jun 18, 2010

@author: miyer
'''

import numpy as np
from collections import defaultdict
#from RNAseq.ordereddict import OrderedDict 

def skip_header(fileh):
    fileh.next()
    for line in fileh:
        fields = line.rstrip().split('\t')
        if fields[0] != '':
            break
    return fields

def read_header(fileh, header_cols):
    header = fileh.next()
    fields = header.rstrip().split('\t')
    col_names = fields[0:header_cols]
    sample_names = fields[header_cols:]
    params = {}
    for line in fileh:
        fields = line.rstrip().split('\t')
        if fields[0] != '':
            # end of params header
            break
        else:
            params[fields[1]] = fields[header_cols:]
    # return the fields from the last line read
    return fields, col_names, sample_names, params

def read_matrix(filename, matrixfilename, header_cols):
    # read file to count number of genes and samples
    fileh = open(filename)
    fields, colheaders, sample_names, params = read_header(fileh, header_cols)
    # count number of genes
    numgenes = sum(1 for line in fileh) + 1  
    # count number of samples
    numsamples = len(fields[header_cols:])
    # close file
    fileh.close()
    # make temp array to store values
    mat = np.memmap(matrixfilename, dtype=np.float, mode="w+", shape=(numgenes, numsamples))    
    # open file again and skip the header
    fileh = open(filename)
    fields = skip_header(fileh)
    # read data into the matrix
    cols = [fields[:header_cols]]
    mat[0,:] = np.array(map(float, fields[header_cols:]))
    i = 1
    for line in fileh:
        fields = line.rstrip().split('\t')
        cols.append((fields[:header_cols]))
        mat[i,:] = np.array(map(float,fields[header_cols:]))
        i += 1
    matobj = type('ExpressionMatrix', (object,), 
                  dict(ngenes=numgenes,
                       nsamples=numsamples,
                       colheaders=colheaders,
                       sample_names=sample_names,
                       sample_params=params,
                       cols=cols,
                       exprs=mat))
    
    return matobj


def create_matobj_dicts(expmat_matobj_genes, expmat_matobj_sample_names):
    # create a dictionary for the gene names
    gene_names_dict=defaultdict(list)
    gene_intervals_dict=defaultdict(list)
    samples_dict = defaultdict()
    for i in range(len(expmat_matobj_genes)):
        gene_name = expmat_matobj_genes[i][0]
        gene_interval = expmat_matobj_genes[i][1]
        gene_names_dict[gene_name].append(i)
        gene_intervals_dict[gene_name].append(gene_interval)
    
    # create a dictionary for the sample names    
    for i in range(len(expmat_matobj_sample_names)):
        sample_name = expmat_matobj_sample_names[i]
        samples_dict[sample_name] = i
    
    return samples_dict, gene_names_dict, gene_intervals_dict



'''
def quantile_normalize(mat, indexes):
    # sort the sample expressions individually and average across genes
    vec = np.average(np.sort(mat[:,indexes], axis=0), axis=1)
    # now reassign the expression back to individual samples
    for index in indexes:
        order = np.argsort(mat[:,index])
        mat[order,index] = vec
        
'''



